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Abstract 

Strong stability preserving (SSP) integrators for initial value ODEs preserve temporal 
monotonicity solution properties in arbitrary norms. All existing SSP methods, including 
implicit methods, either require small step sizes or achieve only first order accuracy. It is 
possible to achieve more relaxed step size restrictions in the discretization of hyperbolic PDEs 
through the use of both upwind- and downwind-biased semi-discretizations. We investigate 
bounds on the maximum SSP step size for methods that include negative coefficients and 
downwind-biased semi-discretizations. We prove that the downwind SSP coefficient for lin- 
ear multistep methods of order greater than one is at most equal to two, while the downwind 
SSP coefficient for explicit Runge-Kutta methods is at most equal to the number of stages 
of the method. In contrast, the maximal downwind SSP coefficient for second order Runge- 
Kutta methods is shown to be unbounded. We present a class of such methods with arbitrarily 
large SSP coefficient and demonstrate that they achieve second order accuracy for large CFL 
number. 



1 Introduction 

1.1 Strong stability preservation 

This work is concerned with numerical methods for the initial value problem 

u'{t) = F{u) u{0) = uo, (1) 

where u E 5?" and F : 5?" — > K". Numerical methods for compute a sequence of solutions 
M^, M^, . . . approximating the true solution at times At, 2At, .... We are interested in numerical 
solutions that satisfy the monotonicity property 

ll""ll < \\u"-'\\ (2) 
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in the case of a one-step method, or more generally 

||m"|| < max{||M"-^||,||M"-2||,...,||M«-*^||| (3) 

for a A: -step method. 

A common approach is to assume that the initial value problem Q is monotone under 
forward Euler integration, subject to some step size restriction: 

\\u+AtF{u)\\ < \\u\\ for all u, for < At < AtpE- (4) 

The term strong stability preserving (SSP) is used to denote any method that gives a solution 
satisfying the monotonicity condition ||3l whenever applied to a initial value problem ([T| 
satisfying the forward Euler condition Til. This can be shown to hold under a step size 
restriction of the form 



At < CAtpE, (5) 

where the factor C, referred to as the SSP coefficient, depends only on the numerical method. 
For a recent review of SSP methods see ||2l. 

The SSP coefficient is, for most known methods, not very large; for a broad class of explicit 
general linear methods it is never greater than the number of stages of the method 16], and 
for many classes of implicit methods it is known or conjectured to be no greater than twice 
the number of stages |[5j[T]|8l. No known methods have C > 2s and exhibit higher than first 
order convergence for large CFL numbers. 



1.2 Downwinding 

The study of SSP methods has been motivated by the numerical solution of hyperbolic con- 
servation laws; in one dimension these take the form 

Ut+f{U), = 0. (6) 

Semi-discretization of the conservation law (|6} leads to the initial value problem Q, where u 
is a finite-dimensional approximation of U and F(m) is an approximation to —f{U)x- 

The solution of the scalar conservation law ||6} has the property that its total variation 
does not increase in time. Hence it is desirable for a numerical discretization to satisfy (|3} 
with respect to the total variation semi-norm; such schemes are said to be total variation 
diminishing (TVD). A common approach to development of TVD methods is to employ a 
semi-discretization for which the TVD property holds under forward Euler integration in 
time, up to some maximal time step size; this is just condition ||4|, where || • || is taken to be 
the total variation semi-norm. 

Such semi-discretizations generally are upwind-biased. By considering corresponding 
downwind-biased semi-discretizations, one arrives at an operator F that approximates f{U)x 
and satisfies 

\\u + AtF{u)\\ < II M II for all Af < AfpE- (7) 
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For example, consider the advection equation 

!if + (ix = (8) 
discretized via the upwind and downwind discretizations: 

"KO = -"^if^ = "KO = -"^1^ = (9) 

Here M/(f) is an approximation to U{iAx,t). In this case the conditions |[7| hold with 
AfpE = Ax, corresponding to the usual CFL condition. 

We say a method is strong stability preserving with downwind SSP coefficient C if the 
method gives a solution satisfying ||3} whenever applied to a system satisfying the forward 
Euler conditions Q and l|7} under the time step restriction 

At < CAipE- (10) 

The idea of using downwinding to achieve the TVD property under larger timesteps for 
explicit time discretizations was originally introduced by Shu and Osher 1131 [141 . These meth- 
ods are frequently referred to as methods with downwind biased discretizations; in the present 
work we will refer to them simply as downwind methods for brevity, with the understanding 
that they incorporate both upwind- and downwind-biased discretizations. Optimal explicit 
downwind Runge-Kutta methods of up to fifth order and ten stages were developed in [V2]. 
Further optimal explicit Runge-Kutta and linear multistep schemes, along with an approach 
to efficient implementation of upwind and downwind WENO discretizations, are given in ||3l- 
A theory of necessary and sufficient conditions for Runge-Kutta methods with downwinding 
to be SSP was developed in [4 1 . 

The principal reason for studying downwind methods is that the downwind SSP coeffi- 
cient C is not as restricted as the SSP coefficient C; for instance, explicit Runge-Kutta methods 
can have order greater than four and C > [12J. In the present work, we consider also im- 
plicit downwind methods. Generally speaking, implicit numerical methods are used in order 
to allow the use of timesteps based solely on accuracy considerations and not on stability. 
As discussed above, implicit SSP methods require step sizes not much larger than those al- 
lowed by explicit methods. The principal aim of the present is to determine whether implicit 
downwind SSP methods allow much larger step sizes. 

1.3 Downwind Methods 

Downwind Runge-Kutta methods take the form 

y. = +Atj2 aijF (f""^ + CjAt,yj'^ + Af ^ S,yF (f""^ + cyAf,yy) , 1 < / < s (11a) 

j=i j=i 

u" = u"-^ +AtJ2 bjF (t""^ + CjAt,yj'^ + Af bjP (t""^ + CyAf,yy) . (lib) 

7=1 7=1 

The method is explicit if each stage depends only on previous stages; i.e., if aij = 0, g,y = for 
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Downwind linear multistep methods take the form 

k-l 

u„ - AthFiun) - AthF{u„) = ^ 0ijUn_i,^j + AtPjF{u„^t:+j) + Atf>jP{Un-k+j)- (12) 

y=0 

The method is explicit if /5j. = j6j. = 0. 

2 Summary of main results 

In this section we present the main results of the paper. The proofs are deferred to later 
sections. 

Our first result concerns explicit downwind Runge-Kutta methods. It is known that the 
SSP coefficient C for a broad class of general linear methods (without downwinding) cannot 
be greater than the number of stages [6J. It turns out that the same bound holds for explicit 
downwind Runge-Kutta methods. 

Theorem 2.1. The downwind SSP coefficient C of any first order accurate explicit downwind Runge- 
Kutta method is no greater than the number of stages of the method. 

It is already known (see lH) that the bound C < 1 holds for explicit linear multistep 
methods. Together with the theorem above, this implies that not too much can be gained by 
using downwinding in explicit methods, at least in an asymptotic sense. Hence we consider 
implicit downwind methods. For implicit linear multistep methods of greater than first order 
accuracy (without downwinding), it is known that the SSP coefficient cannot be greater than 
two. It turns out that the same bound holds for implicit downwind linear multistep methods. 

Theorem 2.2. The downwind SSP coefficient of any second order accurate linear multistep method is 
at most two. 

Given the negative nature of the two theorems above, one might expect that bounds on 
the SSP coefficient hold for the downwind SSP coefficient for all classes of methods. It seems, 
however, that this is not the case; for implicit Rimge-Kutta methods it is conjectured that the 
SSP coefficient can be no greater than 2s, where s is the number of stages. However, it turns 
out that implicit downwind Runge-Kutta methods can have arbitrarily large SSP coefficient. 

Theorem 2.3. For any positive finite r, there exists a two-stage, second-order downwind Runge-Kutta 
method with downwind SSP coefficient equal to r. 



The rest of the paper proceeds as follows. The proof of Theorem 2.1 is presented in 
section |3] The proof of Theorem |2.2| is presented in|4] In Section |5| we prove Theorem 2.3 
by construction and analyze the resulting methods. Some final remarks on these results are 
provided in Section [7| 



3 Explicit Runge-Kutta methods 



In this section, we prove Theorem 2.1 which bounds the downwind SSP coefficient for explicit 
downwind Runge-Kutta methods. SSP methods are most conveniently analyzed by writing 
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them in a certain Shu-Osher form; see |[T6l |7l. The corresponding form for for downwind 
methods is introduced in the following Lemma, which is an easy corollary of Theorems 3.5 
and 3.6 of fT6]. 

Lemma 3.1. The SSP coejficient of the downwind Runge-Kutta method ( (TT) is the largest r > such 
that the method can be written in the form 

^ / At \ / At \ 

y, = u"-^ + E Pi] [yj + y^iyj)) + E P'j [yj + y^iyj)) (i^a) 

s / At \ ^ / At \ 

u" = u"-' + E Ps+i,j [yj + yF{yj)j + E P.+i,j [yj + yf{yj)) (13b) 

with all coefficients non-negative: 

Pij,p,j>0- (14) 

In the remainder of this section, we consider explicit methods of the form | [T3) ; i.e., those 
for which 

Pij = Pij = for all i ^ (15) 

Consider the application of a method satisfying |(T3)-|(T5) to a linear problem with 

F(m) = Lu F{u) = Lw, (16) 

where L, L are fixed matrices, and define z = AtL, z = Af L. Then a straightforward calculation 
shows that the solution can be written as Un = tp{z,z)u„-i where 

i/;(z,z) = X:E'r;7(l + ;)^~Yl + ^) with 7^, > if < r < C. (17) 



7=0 /=o 



Theorem 2.1 follows by considering the conditions for first-order accuracy of the method 
and positivity of the coefficients in the form |[T7|. 



Proof of Theorem 2.1 For a given downwind Runge-Kutta method, consider 1 17 1 with r = C. 
Since F(m) w — F(m), the method approximates the true solution to order p if 

i/;(z,-z) =exp(z)+C(zP+i). 

In terms of the coefficients jji, the conditions for consistency and first order accuracy are: 

EE'r;7 = l, 

7=0 /=0 

ttv (7-20 = C. 
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Since the 7's are positive, we have 



C = EE7;7(y-2/)< . max (y-2/) 
;tS/tS ;6(0,s),/e{0,;) 



□ 



Remark 1. It is possible to show in a similar way that the bound C < s holds even for the 
broad class of general linear methods considered in [6J. 

Remark 2. In |7|, tighter bounds on C were computed for specific classes of downwind Runge- 
Kutta methods. The same approach used there could be applied to find tighter bounds for 
classes of general linear methods. 



4 Implicit linear multistep methods 

Lenferink IITOl showed that the SSP coefficient C is no greater than two for implicit linear 
multistep methods of order greater than one. This bound was shown to hold in an even more 
general sense in |5|. In this section we consider the downwind SSP coefficient C for implicit 
downwind linear multistep methods (T2). Our main result, presented already as Theorem 



2.2 states that the same bound proved by Lenferink also holds for downwind methods. Even 



though this is a more general result than those of ITOl lSl, we provide a simpler proof. 



The method 1 12 1 is accurate to order p if 



k—l k 

E ''if + LiPi - = (0 ^ ^ p) (18) 

;=0 ;=0 

and has downwind SSP coefficient C > r if and only if 

^;,|6;>0 (0<7<fc-l) (19a) 

aj-r{^j + ^j)>0 {0<i<k-l). (19b) 

The following lemma, which appeared in |j6|, facilitates the proof of Theorem 2.2 We say 



a method 1 12 1 of A: steps and order p is optimal if there exists no other method of at most k 



steps and order at least p with larger downwind SSP coefficient C. 

Lemma 4.1. Any optimal downwind linear multistep method | [12) has the property that f^j^j = Ofor 
each j. 

Proof of Lemma [iTj Note that the order conditions | (T8) depend only on the difference /5y — ^j, 
while the inequality constraint | |19a^ can be written as (setting r = C) 

Suppose that an optimal method has /5y > j6y > for 7 e /j c (0, 1, . . . , A: — 1) and f>j > fij > 
for 7 e /2 C (0, 1, . . . ,A: - 1). Then for j e h define ^* = - and f>* = 0; for j e Jz define 
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^* = - fij and /5* = 0. Then the coefficients obtained by replacing /5, ^ with j6* satisfy 
( (19) and (18) with a larger value of C, which is a contradiction. □ 



We now prove our main result on linear multistep methods. Similar to the proof of The- 
orem 2.1 this result follows from combining the order conditions and the positivity of the 



coefficients. 



Proof of Theorem 2.2 It is sufficient to prove that C < 2 for any optimal method. Hence we 



consider an optimal method that is at least second order accurate. Applying Lemma 4.1 



we define fi* = + f>j and cry = sgn(^y - pj) so that f>j - = CTjfi*. Additionally define 
Sj = dj — rfi* . Then the order conditions for order two can be written (taking r = C) 

'■£Sj + Cfi*=l (20a) 

j=o 

Li^i + (Cj + = k - cr^f>l (20b) 
;=0 

'^Y.fSj + {Cf + 2icrj)fi* = k{k - 2cTkfil). (20c) 

Taking -k'^x | [20a) +2k x | [20b) - pOc) gives 

E -{k- + (-C{k - if + 2aj{k - i)) = 0. (21) 

Since all coefficients Sj, /5* are non-negative, at least one of the terms multiplying them must 
be non-negative, for some value of j (else the sum would be negative). Since the terms multi- 
plying Sj are all negative, this implies that 

-C{k-jf + 2aj{k-i) >0 

for some 7; i.e. 

2aj{k-j) 
C < max —7 -w- < 2. 

je[0,k-i] {k-]f 

□ 

Remark 3. Using the order conditions and the SSP conditions, it is straightforward to deter- 
mine (for any particular k, p) the methods with largest C. This was done in [6|. 



5 Implicit Runge-Kutta methods 

It is well known (see, e.g. IIT5l l9ll2l) that (implicit) Runge-Kutta methods of order higher than 
one cannot have C = 00. Furthermore, it is conjectured that they cannot have C > 2s, where s 
is the number of stages IlJ|8l. 
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In im, the following question was posed: do there exist high order implicit downwind 
Runge-Kutta methods with C = oo? Although we have not answered this question directly, 
we have found downwind methods with arbitrarily large C. One class of such methods is the 
two-stage, second-order family 



y2 



r{r - 2) 
At 

yi + 



r 
At 



f(yi) 
f(y2). 



yi 



At 



■4r + 2 



r{r - 2) 



y2 



—Piy2) 



(22a) 

(22b) 
(22c) 



Here we have written the method in the form (13) so that the downwind SSP coefficient is 
apparent: C = r, for any r > 2 + \/2. It can be shown that the method 1 22 1 is A-stable (i.e., 
imconditionally stable when F(m) = —F{u) = Am with 5?(A) < 0). 

Numerical searches indicate that other two-stage, second-order methods with large C exist. 
An obvious drawback of method | [22) is that it is fully implicit. However, a search for two- 
stage, second-order diagonally implicit methods yielded no results. 

For reference, we include here also the coefficients of the method when written in form 



r-2 

r 

r-1 
r 


^ 
^ 


n r2-4r+2 
^ 2r 
n r2-4r+2 
^ 2r 




r-2 1 
2 r 


^ 2r 



(23) 



The left part of (|23j corresponds to the usual Butcher coefficients, and the rightmost part 
displays the additional coefficients Uiphj. In the case F = —F, the method reduces to an 
ordinary Runge-Kutta method; we refer to this method as the underlying method: 



r-2 


r^-2r-2 


r^-ir+2 


r 


2r 


2r 


r-1 


r-2 


r^-ir+2 


r 


2 


2r 




r-2 


4-r 




2 


2 



(24) 



Although the downwind method | [23) will in general behave differently from the underlying 
method | [24) , analysis of the latter may still give useful insight. Straightforward calculation 
shows that the stability function of method | [24| is 



1 + 7^ + W 



r2-4r+2_2' 
2r2 ^ 



(25) 



Further calculation reveals that the method is A-stable since \^{z)\ < 1 for all z with non- 
positive real part. It is clearly not L-stable, but it is nearly so for large r since 

lim ip(z) = — ^ 



is quite small for large r. 
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(a) CFL=0.9 



(b) CFL=8.0 



Figure 1: Square wave advection with the downwind Runge-Kutta method (22 1, compared to 
two standard methods. 



6 Numerical tests 

In this section we conduct a few numerical tests to study the accuracy of the family of meth- 
ods introduced in the last section. We compare the backward Euler method, the implicit 
trapezoidal Runge-Kutta method, and the downwind Runge-Kutta method ( [22) . In all tests 
we take r = 8. 

In the first test we will see that the temporal error for the downwind method includes 
a term related to the spatial discretizations. Consequently, some care must be taken when 
choosing the spatial discretizations. The second and third tests demonstrate that the method 
performs well when appropriate spatial discretizations are chosen. 



6.1 First-order upwind advection 

As a first test we solve the advection equation (|8) using the upwind and downwind differences 
(|9| with Ax = 1/128. We consider the domain < x < 1 with periodic boundary conditions 
and initial condition m(x, 0) = 1 — H(x — 1/2) where H(x) is the Heaviside function. In Figure 
|l]we plot the solutions at t = 1. 

Figure 1(a) shows results obtained with At = 0.9Ax. As expected, backward Euler is more 
dissipative than the second order trapezoidal method. Surprisingly, the (second order accu- 
rate) downwind Runge-Kutta method is more dissipative than even the first-order backward 
Euler method. Figure 1(b) shows results obtained with At = 8.0Ax, the SSP limit for the 
downwind method. In this case, the trapezoidal method, which has SSP coefficient C = 2, 
generates overshoots. The backward Euler method is much more dissipative for this large 
CFL number. The downwind method is more accurate than backward Euler and maintains 
maximum norm monotonicity, but the high level of dissipation exhibited is still disappointing. 

The behavior of the downwind method can be understood through the following analysis. 
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Write the upwind and downwind discretizations as F{u) = Lu,F{u) = —Lu, where L, L are 
the matrices 



/-I 



1 

Ax 



V 



1 



1 \ 



-V 



/ 1 



1 

Ax 



(26) 



V-1 



1/ 



Any downwind Runge-Kutta method applied to this problem results in a recurrence of the 
form (see section 3 above) 



w" = ip{AtL,AtL)u 



n-l 



(27) 



For the method ( 22 1, we find 

, 2 

^p{AtL,AtL) = 

1 + AiL 



- 4r + 2 . 

1 + ^ AtL 

2r 

- 4r + 2 



- 2r - 2 . „ - 4r + 1 . ,2 



2r 



-AtL 



2r2 



-At^LL 



2r 



At (L-L) +C(Af2). 



Since Lm w — Mx/ the first two terms ensure that the method is at least first order accurate. 
However, the term rnvolvrng (L — L) is problematic, since 



/-2 1 



L-L 



1 

Ax 



-2 1 



Vl 



-2/ 



Ax Wv 



(28) 



Thus, if F, F are spatial discretizations with order of accuracy q, then the one-step error for 
the downwind Runge-Kutta method | [22) will contain a diffusive term of 0{rAx^At). In order 
to avoid loss of accuracy, the spatial discretization should ensure that rAx"? ^ Af, so that this 
term is no larger than the 0{Af-) term. 



6.2 WENO advection 

To demonstrate that proper accuracy is acheived when the foregoing condition is satisfied, 
we use a fifth-order WENO interpolation to determine the fluxes. Table [T| shows convergence 
results in the maximum norm for advection of a sine wave (m(x, 0) = sin(27rx) using a CFL 
number of 8. Notice that the downwind method achieves its design order of 2, and gives accu- 
racy similar to the trapezoidal Runge-Kutta method. It should be noted that the downwind 
method is more expensive computationally, since it is not diagonally implicit and requires 
both upwind and downwind operator evaluations. However, the purpose of this test is only 
to confirm the theoretically predicted accuracy. 
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N 


Backward Euler 


2nd order RK 


Downwind RK 


32 


0.728 


0.730 


0.436 


64 


0.603 


0.215 


0.168 


128 


0.452 


0.054 


0.043 


256 


0.292 


0.013 


0.011 



Table 1: Max-norm errors for advection of a sine wave using fifth-order WENO discretization 
and a CFL number of 8. 



6.3 Burgers equation 

In order to provide an initial assessment of the suitability of these methods for application to 
nonlinear hjrperbolic conservation laws, consider Burgers equation: 

Ut+(luA =0 



on the unit interval with initial condition U{x,0) = sin(27Tx) and periodic boundary condi- 
tions. Fifth-order WENO interpolation is again used to determine the fluxes. We consider the 
solution at time t = 0.16, just after a shock has formed. Figure |2] shows a closeup around the 
shock for the exact solution (obtained using characteristics) along with solutions computed 
by each of the methods compared above, using CFL number 6.5 in each case. A grid with 
Ax = 1/512 points is used. As expected, the backward Euler method is much more dissipa- 
tive than the second order methods. Interestingly, the trapezoidal method solution does not 
show oscillations per se, but clearly has an unphysical bump behind the shock. The downwind 
method is by far the most accurate. Again, it should be noted that the downwind method is 
more expensive, since it involves an implicit solve of twice as many equations. 

Figure |3] shows a comparison of two solutions obtained with the downwind RK method 22 



using two different CFL numbers. The solutions are almost indistinguishable, suggesting that 
the spatial error is dominant, as one might expect. Close inspection reveals that the solution 
using the larger CFL number is more accurate. Importantly, it seems that the method does 
not become more dissipative at large CFL numbers. 

We note in passing that solution of the nonlinear system of equations for implicit WENO 
schemes is a significant challenge, especially for large CFL number and in the presence of 
shocks. The CFL number 6.5 used here allowed the use of an easily accessible nonlinear solver, 
namely the newton_krylov and f solve functions of the SciPy package. Efficient solution of 
the nonlinear system for larger CFL numbers is an area for future research. 



7 Discussion 

The new bounds on the downwind SSP coefficient proven here for explicit Runge-Kutta meth- 
ods and implicit linear multistep methods are disappointing, since they indicate that nothing 
can be gained - in an asymptotic sense - by including downwind-biased discretizations for 
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Figure 2: Closeup view of shock in solution of Burgers equation with various methods using a 
CFL number of 6.5. 
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Figure 3: Closeup view of shock in solution of Burgers equation obtained with the downwind 
RK method (|22| with r = 8 and using two different CFL numbers. 
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integrators in these classes. Previous attempts to find methods that have large SSP coefficients 
have been similarly disappointing IITTl lSl. 

The family of second order implicit downwind Runge-Kutta methods we report here is 
therefore quite remarkable. These are the first methods known to have large SSP coefficient 
and to provide higher than first order accuracy in practice for CFL numbers larger than unity. 
They may be useful for problems with hyperbolic components requiring large time steps and 
the avoidance of spurious oscillations. 

In light of this discovery, the following questions are naturally of interest: 

• Do higher (than 2nd) order downwind Runge-Kutta methods exist with large down- 
wind SSP coefficient? 

• Do diagonally implicit downwind Runge-Kutta methods exist with large SSP coeffi- 
cient? 

• Can these methods be implemented efficiently in combination with high order space 
discretizations and large time steps, for problems with shocks? 

Preliminary investigation of the first question indicates the answer is affirmative. Work on all 
three questions is ongoing and will be presented elsewhere. 
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